Aggregate formation in a system of coagulating and fragmenting particles with 

mass-dependent diffusion rates 



R. Rajesh,^'!^ Dibyendu Das,^ Bulbul Chakraborty,-^ and Mustansir Barma^ 

^Department of Physics-Theoretical Physics, University of Oxford, 1 Kehle Road, Oxford 0X1 3NP, UK 
''Martin Fisher School of Physics, Brandeis University, 415 South Street, Waltham, MA 02454, USA 
^Department of Theoretical Physics, Tata Institute of Fundamental Research, Homi Bhabha Road, Mumhai-400005, India 

(Dated: February 1, 2008) 

The effect of introducing a mass dependent diffusion rate ~ m~" in a model of coagulation with 
single-particle break up is studied both analytically and numerically. The model with q = is known 
to undergo a nonequilibrium phase transition as the mass density in the system is varied from a 
phase with an exponential distribution of mass to a phase with a power-law distribution of masses 
in addition to a single infinite aggregate. This transition is shown to be curbed, at finite densities, 
for all a > in any dimension. However, a signature of this transition is seen in finite systems 
in the form of a large aggregate and the finite size scaling implications of this are characterized. 
The exponents characterizing the steady state probability that a randomly chosen site has mass m 
are calculated using scaling arguments. The full probability distribution is obtained within a mean 
field approximation and found to compare well with the results from numerical simulations in one 
dimension. 
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I. INTRODUCTION 

Systems far from equilibrium can undergo phase tran- 
sitions between two types of steady states when the pa- 
rameters of the system are varied. It is important to 
ask about the sensitivity of such nonequihbrium phase 
transitions to changes in the governing dynamics. If the 
transition survives, is the universahty class affected? If 
the transition is lost, does a signature of the lost phase 
remain in any form? 

In this paper, we investigate these questions within a 
lattice model of coagulation and fragmentation in which 
the diffusion constant for a mass m varies as to^" with 
a > 0. For the case in which diffusion is independent of 
the mass {a — 0) and fragmentation involves only chip- 
ping off of unit masses, it is known that there is a phase 
transition from a low-density phase with an exponen- 
tial distribution of masses to a high-density phase with 
a power-law distribution of masses in addition to an infi- 
nite aggregate with a mass proportional to the volume V 
101 . This transition is characterized by a new universality 
class, different from familiar classes such as directed per- 
colation or the parity-conserving class |^ , wetting transi- 
tions, roughening transitions or boundary-driven transi- 
tions 1^ . We will show below that this high-density phase 
is lost as soon as a is nonzero. Remarkably, though, an 
imprint of the infinite aggregate remains in the form of 
a large aggregate that strongly modifies the finite-size 
behavior of the system, and we characterize the scaling 
implications of this. 

Let us summarize the results of earlier related work. 
Enhancement of aggregation moves with increasing mass. 
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corresponding to negative values of a, have been investi- 
gated earlier in the context of coalescing branched poly- 
mers. For a = — 1, using a Smoluchowski approach it 
was shown that the system undergoes a gelation transi- 
tion, i.e., an aggregate that subsumes a finite fraction of 
the total mass forms at finite time. The fragmentation 
move was shown to modify the mass distribution power 
law exponent at the gelation transition B . An off-lattice 
version of the a = case was studied B using Smolu- 
chowski rate equations in the context of aggregation in 
dry environments. In these studies ^, ||, the coagu- 
lating and fragmenting masses represented polymers in a 
solution, undergoing polymerization and depolymeriza- 
tion. In a realistic situation, it may be expected that the 
diffusion of the polymers would depend on their masses. 
For example, in the well known models of polymer mo- 
tion such as the Rouse model or the Zimm model [0, 
the polymer diffusion constant D{m) ^ m^^ and m^^^ 
respectively. This would correspond to a = 1 and 1/2 in 
our model. This provides a further motivation for study- 
ing the model with a mass-dependent diffusion rate. 

Other modifications of the dynamics of the a = 
model that have been studied include changes of the frag- 
mentation rule, the introduction of a spatial bias in the 
dynamics, and the effects of quenched disorder. Intro- 
duction of a mass-dependent fragmentation by allowing 
fractions of masses to break off (as opposed to single- 
particle break up) was studied in @, |[. In this case, it 
could be inferred that the phase transition is curbed in 
all dimensions. Spatial bias was introduced by choosing 
rates such that masses have a preferred direction of mo- 
tion, but with mass-independent hopping rates. In this 
case, it was shown that the phase transition is curbed 
in one dimension In two and higher dimensions, it 
was shown that bias is irrelevant at least as far as the 
existence of a phase transition was concerned. Finally, in 
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a disordered model where fragmentation of masses could 
occur only at fixed sites, it was shown that even in the 
limit of very low disorder, a new mechanism for the for- 
mation of localized infinite aggregates sets in jl^ . 

The remainder of the paper is organized as follows. 
Section II contains the definition of the model, a brief re- 
view of earlier results and a summary of results obtained 
in this paper. Section III contains the analytical proof 
for the nonexistence of a phase with an infinite aggregate 
at large densities for any non-zero value of a. In Sec. 
IV, the exponents associated with the probability dis- 
tribution P(m) are determined using scaling arguments. 
Results of Monte Carlo simulations in one dimension are 
also presented. In Sec. V, the full distribution is ob- 
tained from a mean field approximation and compared 
with the P{m) obtained from numerical simulations. The 
appendix discusses different limiting cases of the problem 
that are solvable exactly. 



II. MODEL AND RESULTS 
A. The model 

The model is defined on a c?-dimensional hyper-cubic 
lattice with periodic boundary conditions. Starting from 
a random distribution of non-negative integer masses at 
each site, the system evolves in time via the following 
microscopic moves : (1) each mass m hops with rate 
D{m) = m~" to one of its nearest neighbor sites chosen 
randomly (2) with rate w, unit mass breaks off from an 
already existing mass and is transferred to a randomly 
chosen neighboring site and (3) following moves (1) and 
(2), the mass at each site adds up. The mass density p 
is a conserved quantity in the model. 

In one dimension, this model can be mapped [0 onto 
other well studied models of nonequilibrium statistical 
mechanics. By interpreting the masses as interparticle 
spacings, the model is mapped onto a one-dimensional 
hard core lattice gas model with competing short and 
long range hops. Correspondingly, the problem may be 
mapped onto a fluctuating interface with competing short 
and long range moves. The limiting case w = oo re- 
duces to the well studied simple exclusion process Q 
or equivalently to a fluctuating interface governed by the 
Edwards- Wilkinson equation 113]. 



B. Previous results for a — 

The case a = was studied by means of a mean field 
approximation [Q, analytical calculations and nu- 
merical simulations in 0, The results are summa- 
rized below. The steady state single-site mass distribu- 
tion P[m) was shown to undergo a phase transition in 
all dimensions. In the p — w plane, there is a critical line 
Pc{w) = VI + w — 1 that separates two types of asymp- 
totic behavior of P(rn). For fixed w, as p is varied across 



the critical line Pc{w), the large m behavior of P{rn) was 
shown to be 



- m /m * 



P(m) 



p < Pc{w), 
P = Pc, (1) 
infinite aggregate, p > Pc{w), 



where by "infinite aggregate" , we mean a cluster that 
contains a finite fraction of the total mass in the system. 
That is, the tail of the mass distribution changes from an 
exponential decay to an algebraic one as p approaches Pc 
from below. As one increases p beyond pc, the asymp- 
totic algebraic part of the critical distribution remains 
unchanged but in addition an infinite aggregate forms. 
All the additional mass in excess of the critical mass 
condenses into this single cluster and does not disturb 
the background critical distribution. The mathematical 
mechanism giving rise to the formation of the infinite ag- 
gregate at the onset of the phase transition was found to 
be very similar to that of the equilibrium Bose-Einstein 
condensation in an ideal Bose gas. 

Finite size effects in the aggregate phase were studied 
in For a system of size V , the probability distribu- 
tion P(m, V) for p> Pc was assumed to have the scaling 
form 



P{m,V) 



^)+l<5[m-(p-p,)y], (2) 



where the exponent x is a crossover exponent, and the S 
function indicates the aggregate part. The exponents x 
and T were shown to be related by the scaling relation 
x(t — 1) = 1. The exponent r was shown to be 5/2 in the 
mean field approximation [Q; further, numerical evidence 
was presented for the exponent being the same in all 
dimensions. 



C. Summary of new results in this paper 

The principal results obtained in this paper are sum- 
marized below. 

(i) It is shown analytically that there is no phase transi- 
tion at finite density for any a > in any dimension. 

(ii) On an infinite lattice with fixed density p, on assum- 
ing a scaling form 



1 / m 
Pm,p = — / - 
\pi' 



(3) 



where f{y) falls exponentially as y ^ oo, it is shown that 
the two exponents are related to each other by the scaling 
relation 

0(2-r') = l. (4) 
The power law exponent t' is shown to be equal to 
2 - ^ for < a < 2, 



for a > 2. 



(5) 
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Equivalently, 

^=(1 for0<a<2, 

^ \ 1 for a > 2. ^ ' 

(iii) In numerical simulations on a finite one dimensional 
lattice, it is seen that an aggregate forms when the total 
mass in the system is increased beyond a certain critical 
value. By analogy with the a = case, we make the 
assumption that P{m) has the scaling form, 

1 / 777 \ 1 

P{m,V)^—g[—)^-5{m-M-M,), (7) 

where Mc is a F dependent critical mass. It is argued 
that t' = T with x being related to r through x(t— 1) = 1 
as in the a = case. The critical mass is shown to scale 
with system size as 

Mc-F2/(2~"\ for a < 2, (8) 

implying that the critical density pc = M^/V diverges 
with the system size. 

(iv) By means of a mean field approximation, we obtain 
the full probability distribution P{m). The scaling form 
Eq. (0) is seen to hold with the exponents as given in 
Eqs.^) and (|). 




m 



FIG. 1: The variation of P(m) with m for a = 1.0 is shown 
for four different values of density p at fixed lattice size V . 
As density is increased, P{m) changes from an exponential 
distribution to a power law distribution. On increasing p 
further, the power law part remains unchanged while the mass 
in excess of a critical density pc condenses into an aggregate. 
The straight line has a power —1.5. The simulation results are 
for a one- dimensional lattice of size 100 and w = 0.1. In the 
inset, the variation of the power law cutoff with system size 
is shown. The simulation results are for a one-dimensional 
lattice with p — 10.0 and w = 0.1. 



III. ARGUMENTS FOR NO PHASE 
TRANSITION AT FINITE DENSITY FOR a > 

On a finite lattice, on increasing the total mass M 
from zero to large values, the following behavior is ob- 
served in numerical simulations. For small values of M, 
P(rn) is seen to have an exponential tail for large mass 
(see Fig. |l|). As M is increased to a critical value Mc, 
P{m) changes to a power law with a cut off at large to. 
As M is increased beyond Mc, an aggregate forms that 
contains all the mass in excess of Mc- The rest of the dis- 
tribution remains identical to the one at Mc- The power 
law part has a lattice size dependent cutoff (see inset of 
Fig. |l|). All these observations are qualitatively similar 
to the a — Q case. A crucial difference is the fact that 
the power law exponent is seen to be less than 2.0 for 
a > 0. This is a puzzle since finite density would im- 
ply that T > 2. In this section, we prove that there is 
in fact no transition at finite densities in the thermody- 
namic limit. The transition seen in finite size simulations 
is explained by the fact that Mc no longer scales as F (as 
in the a = case), but with a power of V greater than 
unity. 

We show that an aggregate with Mc oc V cannot be 
stable at finite densities by assuming the presence of such 
an aggregate and showing that this leads to a contradic- 
tion. In Sec. Ill A, we study the mass profile as a function 



of distance from the aggregate. Based on our observation 
that at distances far from the aggregate the state of the 
system resembles that at the transition point, we obtain 
exact relations that the critical point should satisfy. In 



Sec. [II B, we derive further exact relations by examining 
the two point correlations. In Sec 
the relations obtained from Sees 



IIIC 



[II A 



we show that 
and 



[II B. when 



put together, imply that there can be no phase transition 
at finite densities. 



A. Reference frame fixed to the aggregate 

In the aggregate phase of the a — model, it is known 
that there exists only one large aggregate in steady 
state; if there were more than one, they would collide 
and coalesce into one. This scenario is verified as well 
in numerical simulations for arbitrary a (the area under 
the aggregate part in the mass distribution being equal 
to l/V). Further, in the limit V ^ oo, the aggregate 
becomes immobile for a > because its mass diverges 
with system size. 

Consider a frame of reference that is attached to this 
aggregate. Let TOx and Sx denote the mass and occupa- 
tion probability at a site x with respect to the aggregate. 
Then, by examining the inflow and outflow of mass at 
each site, we obtain 



dt 



= - [wSx + (t- 



'>(l-'5x,o) 



(9) 



with So — 1 and (to^) = J2m=i ^i^)^^ ■ steady 
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state, the time derivative is set to zero. Then, the solu- 
tion of Eq. (|) is 



w for X 7^ 0. 



(10) 



At distances far away from the aggregate, the state of 
the system resembles that at criticality. Taking the limit 
|x| ^ CX3 in Eq. (]To|), we obtain 

{m^-°')^^w{l- Sc) for a > 0. (11) 

This is a relation that the system should satisfy at the 
critical point. 

In the case a = 0, the aggregate is mobile. When 
the aggregate hops, this corresponds to all the other par- 
ticles simultaneously making a hop with respect to the 
aggregate. An analysis, similar to the one carried out for 
a > 0, yields 



2pc — w(l — Sc) for a = 0. 



(12) 



The origin of the factor 2 may be traced to the fact that 
the aggregate is mobile. 



B. Two point correlations 

In this section, we derive further exact relations that 
the system satisfies at the critical point by studying the 
two-point correlations. The analysis is similar to the 
analysis done for the a = case |l^. In the rest of 
the paper, we will work in a coordinate system fixed to 
an arbitrary fixed site. To fix notation, x' will always 
denote one of the 2c? nearest neighbors of x, while Xq will 
denote a neighbor of the origin 0. Let ri{x,x',t) be the 
mass transferred from site x to x' at time t in a time 
interval At. From the definition of the model, it follows 
that 



77(x,x',i) 



with prob. 

with prob. 
otherwise. 



1 At 

2d 



(13) 



To order At, the only nonzero two point correlation in 
the noise is 

(^(xi,x;)2) = ^ [„,l-'^+y,(^i.S,r.^^,o)] ■ (14) 
The mass m^{t) at lattice site x at time t evolves as 



mx(t+A<) =mx(t)-^r;(x,x',t)-H^?7(x',x,t). (15) 



To obtain the two point correlations, we multiply mx(t + 
At) by mo{t -\- At) and take averages over the possible 
stochastic moves and then over the steady state ensem- 
ble of states. Dropping all time derivatives and using 
Eqs. (|l|), (ll) and we obtain 

C„(x) - wD{x) - i- ^ [C„(x') - wD{^)] 



[Ca(0) + WS] ^(5x.O - ^ 51 j ' 



(16) 



where Cq(x) = (rnxWo"") and -D(x) = (mx^mo.o)- The 
homogeneous part of Eq. (|l^ ) is the Laplace equation 
[Cq(x) — wD{-x.)] — 0. With the inhomogeneous part, 
the unique solution is 



Caiji) = w[D{-x) - s] for X 7^0. 



(17) 



Equation (|17| ) is a relation between two point correla- 
tions. A relation between one point functions is obtained 
by summing over all x, the simplifying factor being that 
total mass is conserved. Thus, 



) = wM{l - s) - wsV + ■ 



(18) 



This is an exact relation in all dimensions. 

We are interested in the limit when AI, V ^ oo keeping 
the density p fixed. Taking this limit in Eq. (p^), we 
obtain 



V 



= wp{l-s)-ws, y>l. (19) 



In the exponential phase, (m^~") is finite and hence 
(m}~^)IV ^ as y ^ oo. At the transition point and 
in the aggregate phase (m^^") can at most diverge as 
y^^" (cf. discussions in the later sections of this paper). 
This implies that {m?~°') /V — + as F ^ oo for all finite 
densities and any a > 0. Thus, another exact relation at 
the critical point is obtained: 



= Wpc{l - Sc) - WSc 



Proof of no transition 



(20) 



We combine the results of Sees. Ill A and III B to show 



that there is no transition. The three quantities (m )c, 
Pc and Sc have to simultaneously satisfy two relations, 
namely Eqs. (|ll|) and (|2^). For non-zero values of a, 
this is possible only when cither w — Q and pc = or 
Pc ^ oo. Equations (11) and (20) cannot be satisfied at 



finite nonzero values of pc. This completes the proof that 
there is no transition for a > at finite critical density 

Pc- 

As a check of correctness, pc and Sc can be calculated 
for the a — Q case from Eqs. (^ and ( |20|). W e obtain 
Pc{w) = \/l + 1 and Sc — (w + 2 — 2y/l + w)/w. Not 
surprisingly, this is the result that had been obtained in 
p3[ for the a = case. 

The fact remains that a single large aggregate is seen 
in simulations on a finite lattice when the mass is large 
enough (see Fig. |l]). This observation would be consistent 
with the above result that there is no transition, provided 
the critical density pc seen in simulations diverges with 
V as with /3 > 0. We address this in the next section. 
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IV. SCALING FORMS FOR PROBABILITY 
DISTRIBUTION P{m,p,V) 

A. Large finite densities 

In this subsection, we consider the case when the total 
mass M and the system size V are increased to infinity 
keeping the density p = M/V fixed. In this case, the 
system is always in the exponential phase. We assume 
the following scaling form for the probability distribution: 



lim P(rn, V, p) 



1 



ri fa 



(21) 



where r' and (j) are two unknown exponents. Rigorous 
upper and lower bounds can be placed on r'. Clearly, 
(m) = / dm mP(rn, V, p) should diverge as p when p 
oo. But, the different moments of m vary with p as 



dm m^P(m, p) — dm m^ ^ /< 



^ (22) 

This implies that r' < 2. Also, from Eq. (|9|), (m^-") is 
seen to be finite for all p, in particular for p ^ oo. This 
implies that t' > 2 — a. Also, from the requirement that 
probability distribution sums up to 1, r' necessarily has 
to be greater than 1. These bounds can be summarized 
as 



max(2 - a, 1) < r' < 2. 



(23) 



The two exponents r' and 4> can be expressed in terms 
of one another by an exponent equality. The average 
mass (to) = p. This implies that 



0(2-r') = l. 



(24) 



Thus, there is only one independent exponent, r' is deter- 
mined in Sec. IV B by studying the finite size corrections 



to the probability distribution. 



B. Aggregate formation on large finite lattices 

For a system on a finite lattice, we see (Fig. ^ that 
when the total mass is increased beyond a critical mass 
Mc{V), the probability distribution has a X^-dependent 
cutoff. Any additional mass that is added aggregates to- 
gether to form one massive aggregate. Using this infor- 
mation, we assume the following form for the probability 
distribution: 



Pi-^^V)^±g4^) + ^Sim- 



V 



(M-M,)), (25) 



where M is the total mass in the system. The two expo- 
nents T and X can be expressed in terms of the two other 
exponents r' and (j). We then determine r by scaling 
arguments, thus solving for all the exponents. 



In [Q, it was shown that r and x are related by the 
scaling relation 



x(t-i) = i. 



(26) 



The derivation of this result was based on the fact that 
the number of aggregates is of order unity. The argu- 
ments carry forward to the general a case without any 
modification. We now argue that t' = t from Eqs. ( |2l| ) 
and ( p5| ) . The system feels the presence of the finite size 
when the density dependent cutoff in Eq. ( pl| ) becomes 
of the same order as the lattice size dependent cutoff in 
Eq. (H). That is, when pf ~ V^, or ~ V^^'l'. But, 
is the mean value of the mass in the power law part and 
from Eq. p^^ yx(2-r) rj.^^^^ 



X 



x(2-r) = ^ 



(27) 



Substituting for (j) in terms of r' (see Eq. (p^), we obtain 

r' = T. (28) 

That leaves only one undetermined exponent in terms of 
which all the other exponents can be expressed. 

To determine this exponent, we start with Eq. ( p^ ) at 
the transition point, namely, 

A4(to1-")c-(to2-")^ = wM^{l-Sc)-wScV+wSc. (29) 

Unlike the scaling M/V = p that we used in deriving 
Eq. (^ from Eq. (|l|), we now assume that Mc scales as 
some power of V, namely ^ V'^+^, with /? > 0. From 
Eq. (p5|), we obtain 



/? = x(2-r). 



(30) 



Firstly, by substituting Eq. ( p5| ) in Eq. (p^), it is easy to 
derive that, to leading order in V, {m^~"')c = w(l — Sc)- 
Now, to satisfy Eq. (p9h, there are two cases we have to 
consider : (A) (to^^")c ^ V, or (B) (to^^")c ~ const and 
{m^~°-)c — w[l — Sc) — wScV^^ + . . .. Case (A) requires 
that x(3 — a — r) = 1, which when simplified implies that 
T — 2 — a/2. Case (B) requires that x(3 — a — r) < 
and x(2 — a — r) < — ft which implies that r > 3 — a and 
T > 2 — a/2. For a < 1, these bounds arc in contradiction 
with the rigorous bounds Eq. (psf). Thus for < a < 1, 
only case (A) is viable and hence t — 2 — a/2. For 

1 < a < 2, we have to consider case (B) also. However, 
any solution that arises from choosing case (B) would 
imply a non monotonic dependence of r on a. However, 
we expect that r is a monotonic function of a, and hence 
we discard the solutions arising from case (B). Thus, r — 

2 — a/2. For a > 2, this solution is in contradiction with 



the rigorous lower bound Eq. (23). Therefore, we assume 
that the exponent value is stuck at 1 for all a > 2 (There 
is no contradiction with the above derivation since if the 
distribution were indeed a power law, then the integrals 
would now diverge at the lower cutoff too). This agrees 
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with the exact solution of the a 
in wliich case r = 1. Thus, 



oo case (see Appendix) 



1 for a = 0, 

2 - f for < a < 2, 
1 for a > 2, 



(31) 



where the value for a = is from . Solving for the 

other exponents, we obtain, for < a < 2, 



X 
f3 



2-a' 
a 

2^' 
2 



(32) 
(33) 
(34) 



Now that all the exponents are known, we return to the 
behavior of the scaling function associated with P{m) 
at large finite densities. Numerically, we observe that 
the scaling function fa{x) ^ const as a; ^ for < 
a < 2. For a > 2, we expect fa{x) to go to zero as 
some power of x as a; ^ (see Sec. ^ for numerics). 
This means that, for a < 2, in the limit p — )■ oo, the 
probability distribution is a power law despite the mean 
mass diverging. These observations are consistent with 
the exact solution of the a = oo case (see Appendix) . The 
formation of a power law in the limit of p ^ oo is similar 
to observations in models of aggregation in the presence 
of a constant influx of particles from outside ||l4 , |l^ . In 
these models, despite the mean mass diverging with time, 
P{m) develops into a power law distribution. 

An implication of the exponent r being less than two 
is that the average time scale in the system may become 
very large. The average time scale goes as the average of 
the inverse of the diffusion constant, i.e. (l/_D(m)) = 

(to") ~ TO^"^^ ^ where to* is the mass cutoff ~ p'^ . 
Thus for a > 2/3, it would diverge with m*. On the 
other hand, the inverse of the average diffusion constant 
1/ {D{m)) remain finite, since {D{m)) — (to~") is always 
finite. Thus our model produces a broad distribution of 
time scales with dissonance of average of its inverse, and 
inverse of its average. Such a scenario is reminiscent of 
diffusion in heterogeneous environment which arises in 
super-cooled liquids In the latter system, transla- 

tional diffusion constant averaged over several heteroge- 
neous regions falls out of proportionality with inverse of 
the average time scale. However, the connection of our 
model to the super-cooled liquids should not be taken 
too seriously since while the latter is in equilibrium, our 
model exhibits a nonequilibrium steady state. 



C. Numerical checks 

In this subsection we provide numerical support for the 
assertions in Sec. IV B, from Monte Carlo simulations in 
one dimension. 



Due to finite size effects, it is difficult to 
make an accurate direct measurement of the exponent r 
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FIG. 2; The power law part of P{m) obtained from Monte 
Carlo simulations is shown for three different values of V. The 
simulations are on a one dimensional lattice with w = 1.0, 
p — 15.0 and a — 0.5. The straight line has an exponent 
— 1.75 (see Eq. (|3l|)). In the inset the scaling plot of these 
curves are shown when scaled as in Eq. (Ba) 
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FIG. 3: The variation of pc(V) with V is shown for q = 
0.3 and q = 0.7, where the upper curve has been shifted 
downwards for clarity. The straight lines are best fit with 
power law exponents equal to 0.178 ± 0.004 for a = 0.3 and 
0.539 ± 0.011 for q = 0.7. These values should be compared 
with the analytic results 0.176 .. . for a = 0.3 and 0.538 . . . 
for a — 0.7 (Eq. (B3)). The simulation was done on one- 
dimensional lattices for w — 1.0. 



from Monte Carlo simulations for P{m). However, we 
show that the analytic results for the power law expo- 
nents are consistent with the numerically obtained P(to). 
In Fig. ^, the results from simulations are compared with 
the analytic results for a = 0.5. In the inset, when the 
plots for different V are scaled as in Eq. (p5|), the curves 
lie on top of each other. For a = 1.0, the predicted ex- 
ponent 1.5 also matches very well with simulations (see 
Fig. 0). 

As a second check, we measured Pc(V) as a function 
of y for a = 0.3 and a = 0.7. We adopted the follow- 



7 



ing procedure for measuring Pc{V). We start the sys- 
tem with total mass much greater than the critical mass 
Pc{V)V. The system is allowed to reach the steady state. 
The cluster with the largest mass is identified as the infi- 
nite cluster. PciV) is obtained by measuring the average 
mass in the rest of the system (excluding the infinite ag- 
gregate). In Fig. H, we obtain the exponent f3 from the 
slope of a log- log plot of Pc{V) versus V. There is excel- 
lent agreement with the analytically predicted values. 



V. MEAN FIELD APPROXIMATION 



In Sec. IV, the exponents characterizing the probabil- 
ity distribution P(m) were calculated. These exponents 
were independent of dimension and hence should match 
with the mean field exponents. Also, it was observed Q 
in the a = case that the mean field P{m) matched 
very well with the numerically obtained P{m) for all m. 
In this section, the exponents of the probability distribu- 
tion as well as the full distribution are calculated from a 
mean field analysis. The values of P{m) thus obtained 
are compared with the P{m) for small values of a ob- 
tained from Monte Carlo simulations in one dimension. 
From the mean field analysis, we also calculate P{m) for 
those values of a which are difficult to probe by Monte 
Carlo simulations due to the large times required to reach 
the steady state. 

In the mean field approximation, all correlations are 
ignored by setting all joint probability distributions to 
be the product of single point distribution functions, i.e., 
P(mi,mj) = P{mi)P{mj). Under this approximation, 
the P{m)''s evolve in time as 



dP{m) 
dt 



-P{m){ra + w + s' + ws) + wP{rn + 1) 
P{a)P{m~a) 



}P{m-l) + ^ 



-,m > 0, 



(35) 



dPo 
dt 

where s' 



-s'{l~ s)-ws{l- s)+wPi + s', (36) 

^ij^m~". In the steady state, the time 
derivatives vanish. Multiplying by e"*"" and summing 
m from 1 to oo, and eliminating Pi, we obtain 



Q 



sQ' + ws{l - s){l - e-P) 



ss 



Q' 



w -I- weP -I- wse p 



(37) 



where Q = X^r^Me"^™, Q' = Y."^ P{m)m-°'e-P'^ 
are generating functions. 

The unknown quantities s and s' are determined by 
the two conditions 



(Q')p=0 

dQ\ 
dp) p=^ 



= -P- 



(38) 
(39) 



B 




FIG. 4: P{m) for a = 1.0 obtained from Monte Carlo simula- 
tions (shown in symbols) are compared with the results from 
the mean field analysis (shown as lines). The lattice size is 
V — 400 and w — 1.0 and we have used two densities p = 1.0 
and 2.0. 
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FIG. 5: P(m) obtained from the mean field analysis when 
scaled as in Eq. ( [2l| ) with exponents as in Eqs. (|^ and (||). 
(a) The curves are for a = 1.6, w — 0.1 and for densities 2.38, 
3.81, and 6.69. (b) The curves are for a = 3.0, w = 0.3 and 
for densities 102.74, 165.07, 250.12 and 325.30. The scaling 
function /a (a;) goes to zero for small x. 

Using the series expansion Q = X)n=o("^") ("P)"/""' ^^'^ 
Q' — X]n=o(™" ^-iid comparing terms order 

by order in p, we obtain relations between moments of 
Pirn). From the term in we find 



p{m^ ") = ^^(l — s) — ws. 



(40) 



Interestingly, Eq. ( pO| ) is identical to the exact Eq. ( pj| ) 
in the V oo limit. For a = and a — 1 this yields the 
two results s = -^^^r^ pi and s = /'"^ . 

Comparing the terms proportional to , we obtain the 
relation 



p{m 



ws 
P 



pw{l + s). 



(41) 



Thus, for large p, {m ) ^ p {m ) . This provides us 
with a method for deriving the exponents from the mean 
field equations. Assuming the scaling form Eq. (Eh) and 
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using the exponent identity Eq. (|2J) , there remain one in- 
dependent exponent to calculate. Using the scaling form 
in p, we obtain 



(/.(3 - r') = 2 + 



a), 



(42) 
the 



which immediately yields t' = 2 — a/2 and (p ~ 2/q;, 
same as in Eqs. ( pl|) and (|3^). 

We now calculate numerically the full P{m) from the 
mean field Eqs. (|3^) and (^6|). If s and s' are known, 
the full distribution P{rn) is known. We use this fact to 
determine the full distribution numerically by the follow- 
ing procedure. We fix s and s' at a certain initial value 
and calculate the resulting P{m) and check the consis- 
tency condition s' = ^ P{m)/m°'. We tune s' to satisfy 
the above condition to an accuracy of 10~^, to determine 
P(m)'s, and thus the density p. 

Using the above numerical method, P{m) was calcu- 
lated for various a's. In Fig. ^, we compare some of these 
mean field results with P{m) obtained using Monte-Carlo 
simulation, for large V . The agreement is excellent, sug- 
gesting that the mean field results are a very good ap- 
proximation to the actual answer. 

We now use the mean field results to probe P{m) for 
values of a that cannot be studied easily by Monte Carlo 
simulations. In Fig. ^, we show the scaling plots for a ~ 

the small x 
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1.6 and a = 3.0. As mentioned in Sec. 
behavior of the scaling function has a different behavior 
for a < 2 and a > 2. In the former case f{x) ~ const, 
while in the latter case f{x) ~ when a; — > 0. 



If the fragmentation move was to be considered on its 
own, then a fluctuation of order I would be dissipated 
in time of the order P. This exponent is known ex- 
actly because of the exact analogy in one-dimension 
between an only-fragmentation model and the Edwards- 
Wilkinson interface ||l^. For a = 0, the two processes 
are of similar strength and hence there is the possibility 
of a transition. But for a > 0, the fragmentation process 
always dominates and hence there is no aggregate phase. 
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APPENDIX A: EXACTLY SOLVABLE LIMITS 

In this appendix, we discuss the limiting cases of the 
model for which the full probability distribution P(m) 
can be calculated. 



VI. SUMMARY AND CONCLUSIONS 

In summary, we have studied the steady state of a 
system of aggregating and fragmenting particles, with a 
mass dependent diffusion rate D{m) ~ m~" with a > 0. 
We showed analytically that the nonequilibrium phase 
transition which is known to exist for a = 0, vanishes 
when a > 0. This is in agreement with the results of 
extensive numerical simulations, through which we ex- 
plored the dependence both on system size and total 
mass. Although no true infinite aggregate forms in the 
thermodynamic limit, its imprint at high densities re- 
mains in finite-sized systems in the form of an aggregate. 
Further, for the single site mass distribution function, we 
obtained the exact scaling exponents associated with its 
dependence on the mass, the density and the system size. 

Our results give more credibility to the intuitive argu- 
ments presented in as to the circumstances in which 
one should expect to see a nonequilibrium phase with 
an infinite aggregate, as occurs in the a = case. We 
reproduce the argument here. In the model under con- 
sideration, there are two competing processes: while the 
diffusion move creates larger and larger masses by coag- 
ulation, the fragmentation move tends to create smaller 
masses, as well as to inhibit the formation of large masses. 
If the diffusion move was to be considered by itself, then 
a cluster of size / would be created in time of order 



In the limit a — > oo, the rate of diffusion becomes equal 
to zero for all masses m > 2. The model then reduces 
to a zero range process in which with rate w unit 
mass can break off from masses m > 2, while the unit 
mass can hop to a neighboring site with rate 1 + w. It 
is then straightforward to verify that the steady state 
probability distribution has a product form, i.e., 

P{...,mi,m2,...) = \[P{mi), (Al) 



with 



P{m) ^ \ c{i+w) 



C7™ "1 > 1, 



TO = 0. 



(A2) 



The constants c and 7 are fixed by the two constraints 
Pi'fTT') = 1 a-iid J2m ™P(™) = P- Solving for c and 
7, we obtain 
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w{l — s) 
1 + w ' 
s{w + 1) 



w + s 

with the site occupation probability s being equal to 



(A3) 
(A4) 



y^W^l + p)^+4pW-w{l + p) ^^^^ 
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In the limit w — > oo, s has the correct hmit p/(l + p) (see 
Eq. (P)). 

We would be interested in the form of P{m, p) when 
p — > oo. Expanding s in terms of we obtain, 



l + wl 1 , 

5 = 1 

w p 



In this limit 



1 



Pim,p) w -e^"/", 
P 



oo. 



Thus 



P{m,p) ^ —foc{—), m,p~^(X), 
m p 



(A6) 



(A7) 



(A8) 



wher e th e scaling function foa{x) ~ x when a; ^ 0. From 
Eq. (AS), we see that r = 1 for a oo. 



In the limit w = oo, the model reduces to a zero range 
process |17 . As in the a = oo case, the steady state prob- 



ability distribution has a product form as in Eq. (jAl|). 
P{m) for this limiting case was worked out in ||l|. For 
the sake of completeness, we reproduce the final result, 



P(m) = 



1 



TO > 0. 



(A9) 



lu = 



In this limit, masses diffuse and coagulate on contact. 
Clearly, the steady state is one in which the entire mass 
is clumped together into one aggregate. For the a > 
problem, this is the only limit in which an aggregate 
forms which holds a finite fraction (here unity) of the 
total mass at finite density. 
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